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Abstract 

In all kinds of engineering problems, and in particular in methods for 
computational fluid dynamics based on regular grids, local grid refinement 
is of crucial importance. 

To save on computational expense, many applications require to re¬ 
solve a wide range of scales present in a numerical simulation by locally 
adding more mesh points. In general, the need for a higher (or a lower) 
resolution is not known a priori, and it is therefore difficult to locate areas 
for which local grid refinement is required. In this paper, we propose a 
novel algorithm for the lattice Boltzmann method, based on physical con¬ 
cepts, to automatically construct a pattern of local refinement. We apply 
the idea to the two-dimensional lid-driven cavity and show that the au¬ 
tomatically refined grid can lead to results of equal quality with less grid 
points, thus sparing computational resources and time. The proposed 
automatic grid refinement strategy has been implemented in the parallel 
open-source library Palabos. 


1 Introduction 


The lattice Boltzmann method (LBM) is a numerical method which is widely 
used in the field of computational fluid dynamics. It has proven its importance 
with respect to other traditional numerical methods, as it has recently been 


used in many simulations of engineering interest (see Aidun and Clausen 2010 
Parmigiani et al. [2011 , Malaspinas and Sagaut [2012] among others). 
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In order to overcome the limitations of the method which, in its most 
straightforward formulation, is restricted to uniform grids, the community has 


proposed several approaches to block-wise grid refinement (see Filippova and 
Hanel 1998 , Dupuis and Chopard 2003 , Chen et al. 2006| ). The goal when 


using such methods is to increase the precision of the results through a local 
adaptation of the mesh resolution, with an acceptable computational effort. 

In this paper we treat, independently of the underlying grid refinement al¬ 
gorithm, the question of how a grid refinement pattern is to be devised most 
efficiently. We must note that, while this point remains as of now open in the 
lattice Boltzmann community, this operation is of crucial importance, as all the 
potential gains of non-uniform grids are tightly coupled with their choice of 
arrangement in space. 

In traditional numerical methods such as the finite differences, there exist 
frameworks for a posteriori error estimation. This allows to locally mark areas 
where the biggest local errors are found. The interested reader is invited to read 
Skeel 1986 for a list of error estimation techniques. 


In the case of the LBM, this kind of error estimation has not been stud¬ 
ied in depth. We can however cite Han et al. [20061 where the authors adapt 
the Richardson extrapolation method (commonly used in finite difference) to 
the LBM. They find that the error for two grids with different resolutions is 
proportional to the off-equilibrium part of the distribution function. 

Another approach is proposed in Crouse et al. 2003 . The authors use the 
so-called sensors <p s to provide local error estimates. In particular, in a grid 
with spatial discretization Sx , they work with the divergence sensor defined as 


= Sx 3/2 \V ■ 


( 1 ) 


which is compared to an empirical value in order to generate a non-uniform 
grid for their simulation. 

In this work, we propose a novel method to find zones that might require 
a local finer resolution, based on an anticipation of how the ratio between off 
equilibrium and equilibrium parts of the velocity distribution functions con¬ 
verges towards the Knudsen number (Kn). In the convective rescaling, the 
Knudsen number is simply a constant, independent of the level of grid refine¬ 
ment. The novelties of our approach reside in the fact that we do not have to 
perform several simulations to estimate the local errors and that it also provides 
an estimation of the degree of refinement needed. 

The structure of this document is as follows. We start by briefly introduc¬ 
ing the LBM in order to define all the theoretical concepts and notations in 
Section [2] Then, in Section [3| we present the generic concepts of our criterion 
for automatic refinement. A numerical validation for the 2d lid-driven cavity is 
provided in Section [4j and Section [5] provides a conclusion and a discussion of 
future work. 
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2 The lattice Boltzmann method 


The lattice Boltzmann method (LBM) is by now a well-known numerical method 
for computational fluid dynamics. In this section we have chosen to present only 
the basic concepts. For further details the reader is referred to |Chopard and| 
Droz [2005] , Sued [2001] , Wolf-Gladrow [2000] . 

We start by presenting the Boltzmann equation (BE). This equation de¬ 
scribes the time evolution of large numbers of particles in a region of the space 
x £ M 3 with a given velocity £ £ R 3 , which are represented by the particle mass 
distribution function /(:c,£,t). The Boltzmann equation for a gas without ex¬ 


ternal force with the BGK (for Bhatnagar, Gross, Krook, see Blratnagar et al. 
1954 ) approximation with relaxation time r is 


(d t +£ • V*)/ = - 1 (f{x,£,t) - / eq (cc,|,t)), 


( 2 ) 


where / eq is the equilibrium distribution function, given by the Maxwell-Boltzmann 
distribution. 

Following the ideas presented in Shan et al. 2006] , one can discretize the 
velocity space to a set of q velocities. For example, in 2D, a well-known set of 
discrete velocities is given by the D2Q9 model (see Fig. [l]) 


{&}!=o = 0), (-1,1), (-1,0), (-1, -1), (0, -1), 

( 1 ,- 1 ), ( 1 , 0 ), ( 1 , 1 ), ( 0 , 1 )}. 


(3) 



Figure 1: The D2Q9 lattice with the vectors representing the microscopic ve¬ 
locity set . A rest velocity £ 0 = (0,0) is added to this set. 


The BE after velocity discretization becomes 

{dt + & ■ V*)fi(x,t) = -~(fi{x,t) - f- q {x,t)), (4) 

T 

where we have used the notation fi(x,t) = 

Next, this equation is integrated along characteristics with the trapezoidal 
rule. One obtains the following implicit equation 

fi(x + 6t£i,t + 6t) - fi{x,t) = 

- ^{fi( x + S t€i,t + S t) + /»(*, t) 

-fPtx + St^t + S^-fP^t)}. (5) 
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This equation can be converted in an explicit scheme by using the change of 
variable given by 


fi = fi + - /D 

2 t + St 
T= 2 St ’ 


( 6 ) 

(7) 


Finally, we obtain the LBM BGK equation 

fi(x + 5t£i,t + St) - fi(x,t) = ~[fi{x,t) - fP(x,t)}. (8) 

T 

We note that Eq. ([8]) is the kinetic equation that we simulate in the LBM. From 
now on the bars are omitted for simplicity unless explicitly stated otherwise. 
Without loss of generality we also set St = 1. 

The discrete equilibrium distribution function is expressed by the truncated 
Maxwellian equilibrium 


f! q = WiP ^1 + + ^4 Qi : ) ( 9 ) 

where p is the density, u is the macroscopic velocity field, Qi = — c 2 s I, c s 

and Wi the lattice speed of sound and the lattice weights respectively, which are 
given for the D2Q9 lattice by 

c 2 s = 1/3, w 0 = 4/9, IC2,4,6,8 = 1/9, 101,3,5,7 = 1/36. (10) 


The density and the velocity fields are computed by the distribution function 
through the relations 


9 -l q-1 

5> = £/r. 

(11) 

2 — 0 2—0 

q-1 q-1 

(12) 


2 — 0 2—0 


For implementation purposes, a time-step is decomposed into two parts that 
are applied successively on the whole computational domain. The two-steps are 
called the “collide-and-stream” operation. 

1. The collision, which modifies locally the value of the populations according 
to 

fi ut ( x ,t) = M x ,t) - \ ( M x ,t ) - yrW)) • ( 13 ) 

2. The streaming, which moves the populations to their neighbors according 
to their microscopic velocity 

f t (x + £i,t + 1) = f° ut (x,t). (14) 
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If we perform a multi-scale Chapman-Enskog (CE) expansion (see Chapman 


and Cowling 1960 , Chopard and Droz 2005 for more details) one can show that 


the LBM BGK scheme is asymptotically equivalent to the weakly compressible 
Navier-Stokes equations 


dtp + V • {pu) = 0, (15) 

d t u + (u-V)u = --Vp + 2uV ■ (S), (16) 

P 

with p being the pressure, S the strain tensor and v the kinematic viscosity 
defined by 

P = c 2 s p (17) 

S = ^(Vu + (Vu) T ), (18) 

u = c 2 s (t- 1/2). (19) 


The CE expansion is done under the assumption that /) is given by a small 
perturbation of the equilibrium distribution 


fi = fr + efl 1] +0(e 2 ), 


( 20 ) 


where e <C 1 can be identified with the Knudsen number (see Huang 1 987] ). 
This can be seen by replacing the CE Ansatz (Eq. (201) in Eq. (|4j), neglecting 
the time dependence of /, and keeping only the lowest orders on both sides of 
the equation one remains with 


& • v/, 

Cs ,(0) 

L Ji 


(o) 


1 f (i) 

T h ’ 


IfU) 

T h > 


(i) 


( 0 ) 


TCs 

~L~ 


= 1 = K„, 


( 21 ) 


with Kn the Knudsen number, A = c s t the mean free path and L the char¬ 
acteristic length of the system. In the second line we used the fact that the 
characteristic value of the microscopic velocity is given by the speed of sound, 
and that since we are interested only on the variation of macroscopic quantities, 
the gradient must scale like the characteristic size of the domain. The fi can 
also be formally decomposed into 


fi = /r + 


f ne Q 


( 22 ) 


where / neq is the off-equilibrium part of the particle distribution function. Ac¬ 
cording to the assumptions of the CE expansion /” cq can be approximated to 
/W when one neglects 0(e 2 ) terms 

rq-e/W. (23) 
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( 24 ) 


During the CE expansion, one finds that is given by 


E/ ‘ = 2? Qi:Il(1) ’ 


where the tensor ifi 1 ) = JA £i£i e fi^ is related to the strain rate tensor S 


through the relation 


IlW = -2 c 2 sP tS. 


(25) 


3 Grid refinement criterion algorithm 

In this section we present the grid refinement criterion, independently of the 
actual algorithm used for the grid refinement. All relevant details about the 
grid refinement approach used in the validation section of this paper can be 


found in Lagrava et al. 2012 


Our grid refinement criterion is based on the fact that the off-equilibrium 
and equilibrium parts of the distribution function are related by the Knudsen 
number through the relation 


yneq _ jeq Knj 


as shown in Eq. (21). 


Let us now consider a grid G c with spatial spacing 8x c . We also define a finer 
grid Gf with spatial step 5xf , defined by the relationship 8x c = Sxf/n, where n 
is a positive integer. Let’s assume that the same simulation is executed with both 
grids, and that the Reynolds number Re is the same in both cases. We recall 
that Re = UfoN/i/n,, where uib, r'lb, and N are respectively the characteristic 
velocity, viscosity and length scale of the system in lattice units. We also set 
■uib to be constant in G c and Gf (convective scaling), a fact which can be 
expressed through the relationship St c /Stf = 5x c /5xf. This also implies that 
both simulations have the same Mach number 


Ma = uib/c s . 


(26) 


Finally, we note that the Knudsen number, which is defined as 


Kn = X/L, = c s r /L 
ji 


u ib ct s t 
c s Lnib 


Mlb V 
Cs Luib 


= Ma/Re, 


(27) 


is independent of the resolution. As a side remark, this statement would not 
hold in case of diffusive scaling, i.e. when the parameters between two grids 
scale like St ~ Sx 2 . 

The general idea of our refinement criterion is that the quality of our results 
are given by the measured Knudsen number (see Eq. ( |21[ )) when compared with 
the “theoretical” Knudsen number of Eq. (27). To formalize the automatic grid 


refinement process, we propose the following algorithm: 
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1 . 

2 . 

3. 

4. 

5. 


Choose a unique spatial resolution 5x for the computational domain R. 
Divide the simulation domain R in m-by-n sub-domains Rij- 
Perform the simulation over R. 

Compute the Kn number from the relation Kn = Ma/Re. 

For each sub-domain Rij of R compute the mean value of the quantity 


1 1 q ^ 1 

Ci 'i = q |j/7T| 

q 1n ^ xeRij i=o 


ir 

fP 


(28) 


where \R,i..j is the number of grid points inside the Rij region. 


At this point one wants to determine if each region R,,j has a grid has a sufficient 
resolution. To do so we propose to compute the refinement factor R.[ ■ defined 
as 

round (log 2 ^)V (29) 

where roundQ denotes that the result of the logarithm is rounded to the closest 
integer. This quantity determines the ratio between Sx on our test grid R and 
the target resolution on a subdomain Rij one would like to obtain in order for 
its grid to be resolved enough. For example if R{j = 0 then no refinement is 
needed on the subdomain R,j . but if R.j :) = 2 then the grid in Ri j must be 
made four times finer. 


4 Numerical results 


We have chosen to test our algorithm for a steady state problem, as a dynamic, 
time-dependent adaptation of the grid is not the scope of this paper. The test 
case is the well known 2D lid-driven cavity, for which reference benchmark values 


are provided in Ghia et al. 1982 


The numerical setup and structure of the solution are depicted in Fig. [2] 
We work on a square bounded domain with no-slip walls, except for the top lid, 
which is subject to a constant right-directed velocity. In our particular case, the 
Reynolds number is fixed to Re = 100. The discrete time step of the simulation 
is pinned down by setting the velocity of the top lid in lattice units to uib = 0.01. 
In this numerical example, the Knudsen number is Kn = uib/(c s Re) = 0.00017 
(see Eqs. (26) and ©)• 

We divide the domain in several regions R,_j (the example used here with 
five lines and five columns can be found in Fig. [3j|. We perform the measure 
of the quantity for a uniform resolution with N = 15,30,60,120, 240. The 
values of Cij /Kn for N = 30 are depicted in Fig. [4] The required refinement 
factor is then computed by taking the binary logarithm of C 7 ;. ? /Kn. 
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Figure 2: Set-up for the cavity 2D example and norm of the non-dimensional 
velocity at steady state. 
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Figure 3: Division of the simulation domain in several regions. 


Fig E compares the value of the Cjj coefficient to the Knudsen number at 
different levels of resolution. A white color means that the region R l ;j must be 
refined {C,^ > Kn), and gray means that the resolution is sufficient. It can be 
seen that at N = 15, the resolution is insufficient in the whole domain. At larger 
values of N, we observe that the white regions decrease until only the corners 
are left. It should be pointed out that the predictions of our algorithm seem 
reasonable, as the velocity gradients are obviously largest close to the top lid. 
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Figure 4: Fixing N = 30, i?T for each region Ri.j. 


Furthermore, it makes sense to require a higher grid resolution in the corners, 
in which the velocity imposed by the boundary condition is discontinuous. 




N = 60 TV = 120 





Figure 5: Areas in a 5-by-5 grid that require grid refinement, depending on the 
level of grid resolution TV. White color means that further refinement is need 
(i.e. Cij > Kn), while gray stands for the fact that the resolution is sufficient. 

Fig.[6]shows the grid dependence of the largest, the smallest, and the average 
value of the ratio between Cjj and Kn. It can be seen that this quantity obeys 
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a power law with exponent 1. In other words, one can take it for granted that 
Cij is at most reduced by a factor two when the resolution is doubled. 



Figure 6: Convergence of several C it j values divided by Kn. 

Using the results presented in Fig. [4j we generated a non-uniform mesh 
which is described in Fig. [7] The level zero is corresponding to N 0 = 30, level 

Top lid 


2 

2 

2 

2 

2 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

0 

1 

1 

1 

0 

0 

0 

0 

0 

0 


Figure 7: Level of refinement of each block on our non uniform grid. 


one to Ni = 60, and the level two to N 2 = 120. In order to assess the quality 
of our automatic refinement technique, we generated a “naive” mesh (shown in 
Fig. [8]) . We also compared the results obtained with the two non-uniform grids 
with the results of a uniform grid with N = 120, and with reference results of 
Gliia et al. 1982 . In total, the number of points used in the “automatically” 


generated grid is 5436, while there are 6400 points in the naively generated one, 
and 14400 points in the uniform one. The economy in grid points of the naive 
mesh if of 44.4% as compared to the uniform case, and the economy of the 
automatic approach is of 15% as compared to the naive one. 

We have computed the root mean square (RMS) for the centerline ^-component 
of the velocity with the results of Ghia et al. |l982| at a Reynolds of Re = 100. 
The three results for the centerline velocity are depicted in Fig. [9] We have found 
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RMS values which are the same as the RMS value obtained with a uniform grid 
at N = 120, which is given by 

RMS 120 = 0.00208, RMS gr = 0.00213, RMS naive = 0.00219, (30) 

Here, RMS 120 , RMS gr , RMS na i ve are the RMS values of the uniform grid, the 
automatic, and the naive grid refinement strategy respectively. We therefore 
conclude that our criterion for grid-refinement allows for results of good quality 
while requiring less computational power than the naive or uniform grids. 



Figure 8: Three-level refinement, obtained by approximation, to achieve a good 


value of RMS as compared to Ghia et al. 1982 



Figure 9: Vertical centerline u x velocity for reference solution, uniform and 
non-uniform grids. 

Finally, the generation of the grid of Fig. [7] from Fig. [3] is not completely 
straightforward at first glance, since the resolution levels do not match between 
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the two figures. Nevertheless one must keep in mind that for our grid-refinement 
algorithm to be valid, on adjacent subdomains only a factor two is allowed for 
the ratio of resolutions. For example if Rq 0 = N then Rq 0 = N ± 1 can 
actually be implemented. Therefore the inconsistency that seems to appear at 
first between Figs. [7] and [4] is rather due due to this consistency constrain of 
resolution level ratio between different resolution regions. 


5 Conclusions and future work 

In this work we presented a novel algorithm, based on physical arguments, 
which allows a detection of regions where local grid refinement is needed for 
the lattice Boltzmann method. The detection is based on the evaluation of the 
local Knudsen number which is computed from the density distribution func¬ 
tion. This method has the advantage on relying only on local information and is 
therefore very efficient in terms of parallelism. Our method is validated on the 
two dimensional lid driven cavity and shows that without any a priori knowl¬ 
edge of the flow a good mesh is proposed and the results of the simulation are 
matching those of meshes (homogeneous and locally refined) with significantly 
more grid points. This technique is can be straightforwardly applied in three 
dimensions. Since no a priori knowledge of the numerical setup it could also be 
used for adaptive grid refinement techniques in time dependent problems. 
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